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molecular dynamics simulations and bond valence analysis 
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Based on molecular dynamics simulations of a lithium metasilicate glass we study the potential of 
bond valence sum calculations to identify sites and diffusion pathways of mobile Li ions in a glassy 
silicate network. We find that the bond valence method is not well suitable to locate the sites, but 
allows one to estimate the number of sites. Spatial regions of the glass determined as accessible for 
the Li ions by the bond valence method can capture up to 90% of the diffusion path. These regions 
however entail a significant fraction that does not belong to the diffusion path. Because of this low 
specificity, care must be taken to determine the diffusive motion of particles in amorphous systems 
based on the bond valence method. The best identification of the diffusion path is achieved by using 
a modified valence mismatch in the BV analysis that takes into account that a Li ion favors equal 
partial valences to the neighboring oxygen ions. Using this modified valence mismatch it is possible 
to replace hard geometric constraints formerly applied in the BV method. Further investigations 
are necessary to better understand the relation between the complex structure of the host network 
and the ionic diffusion paths. 

PACS numbers: 66.30.Dn,66.30.Hs 



I. INTRODUCTION 

On a coarse grained time scale, ionic conduction in 
glasses is commonly described by a thermally activated 
jump motion of the mobile ions in an irregular host 
framework formed by a network former (Si02, B2O3,...). 
Regions of long residence times in the network are associ- 
ated with sites, and rare events, where ions move between 
neighboring sites, are associated with jumps. The prob- 
lem how to identify and characterize sites and diffusion 
paths in such non-crystalline materials poses a challeng- 
ing task of current research activities. Its solution can 
be expected to play a key role in approaching a more 
quantitative understanding of ion transport properties in 
glasses; 1 Methods developed in this field should be useful 
also for the description of thermally activated transport 
processes in other disordered systems. 

In recent molecular dynamics simulations of modified 
network glassesS^^iL^ sites were identified based on 
the number density p(x) of the mobile ions. In some 
of these works, the sites constitute those spatial regions, 
where p(x) exceeds a certain threshold p*. The threshold 
value was specified by requiring the number of these re- 
gions to be maximal (for details of the procedure and the 
question which of these regions are identified with sites, 
see ref.0and the analysis in Sec. lIIIf) . The sites form dur- 
ing the cooling process in the highly viscous mehpi and 
are stable on the time scale of the main structural rear- 
rangements of the network, which, below the glass tran- 
sition temperature, is much larger than the time scale 
needed for the mobile ions to reach the long-time diffu- 
sion regime. For the systems investigatedj2iiMiiii§iL& the 
number of sites exceeds the number of ions by 6-12% only, 



suggesting that a vacancy mediated hopping dynamics 10 
should be considered in a coarse grained descriptioniiiS 
Due to the interaction effects the ionic motion can be 
strongly cooperative in such situations, as reported by 
one of the authors*^ and also in ref. [l3|. 

Although the method of determining the local number 
density in MD simulations was by following the trajec- 
tories of the mobile ions and by registering their occu- 
pation times in small cells (over sufficiently long time 
periods), it should be noted that p(x) can be viewed as 
an equilibrium quantity with respect to configurations 
belonging to the metastable attraction basin of the free 
energy reached after the cooling process. This means, 
that in the absence of aging effects, which involve tran- 
sitions between different basins, p(x) can in principle be 
determined independently of the ion dynamics. 

Experimentally, a measurement of the local number 
density of mobile ions with respect to the host struc- 
ture of the network forming ions is currently not pos- 
sible. Unlike in crystals, it is already difficult to ob- 
tain a reliable representation of the framework struc- 
ture. The most successful procedure to date is to use 
x-ray and neutron diffraction data as input for a re- 
verse Monte Carlo modeling^ Structures obtained in 
this way are in agreement with the static structure fac- 
tors and they fulfill certain constraints imposed by chem- 
ical requirements. Recently this structural modeling has 
been combined with bond valence (BV) sum calcula- 
tions in order to identify diffusion paths for mobile ions 
in network glasses i 15 i 16 ' 17,18 However, while for various 
crystalline systems it could be demonstrated that both 
sites and pathways predicted by the BV method match 
the results from detailed anharmonic crystal structure 
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analysesji2i22i2L2S the potential of the BV method for 
predicting diffusion pathways or sites in amorphous sys- 
tems needs to be clarified. 

Since ionic transport pathways in glasses cannot be di- 
rectly inferred from experiment, we test the BV method 
in this work by MD simulations. The interaction poten- 
tials used in our MD simulations were derived from ab 
initio Hartree-Fock self-consistent calculations and were 
checked to keep the crystal structure stable under con- 
stant pressure conditions^ This is a severe test of the 
quality of the potential parameters. Many properties 
of ion conducting glasses, as e.g. structural information 
obtained from scattering experiments, vibrational spec- 
tra, values for ionic diffusion and various thermodynamic 
properties are successfully reproduced by the MD model. 
Using the MD simulations as a reference to test the BV 
method, we nevertheless have to bear in mind that the 
underlying Hamiltonian for the interactions between the 
ions and the way of vitrifying in the MD simulations (e.g. 
the high cooling rates) may not give a fully accurate rep- 
resentation of all properties of the real experimental sys- 
tem. In any case, even if the representation with respect 
to the sites and conduction pathways is not perfect, one 
should expect that the BV method is still applicable, 
since the MD model itself presents a valid physical sys- 
tem. For the model system, however, the optimal param- 
eters used in the BV analysis could possibly be different 
from those applied to the corresponding real glass (cf. 
Sec.lTvL 

Our results show that the BV paths entail the sites of 
the mobile ions in the MD simulations. The BV method 
yields a reasonable estimate of the number of sites for the 
mobile ions, but is not well suitable for locating the sites 
within the pathways. In order to evaluate the potential of 
the BV method with respect to diffusion pathways, the 
latter have to be specified from MD simulations. Such 
specification can be done based on p(x), by using the per- 
colation threshold p pcrc for attributing the regions with 
p(x) > pp 0rc to the diffusion path. We will show that the 
BV path reaches a sensitivity up to 56%, i.e. it covers up 
to 56% of this diffusion path, and, conversely, the speci- 
ficity lie in the range 50-60%, i.e. 50-60% of the BV path 
belong to the diffusion path. A variant of the BV method 
recently developed by one of the authors^ allows us, by 
introducing a penalty function for unfavorable asymmet- 
ric bonding situations, to achieve a sensitivity up to 90% 
and a specificity between 30-60%. 

The paper is organized as follows. After a description 
of the MD simulations in Sec. [HI we perform an identifi- 
cation of the ion sites in Sec. IIIII following the procedure 
suggested in ref. 0- We include a discussion of the role 
of the grid spacing used in this procedure and suggest 
an alternative criterion to finally assign the regions with 
p(x) > p+ to sites. The diffusion path is determined 
from a percolation analysis. In Sec. IIVI we perform a 
BV analysis based on time-averaged network configura- 
tions obtained from the MD simulations and evaluate in 
Sec. |V] the potential of the BV method for identifying 



TABLE I: Potential parameters for the MD simulations (cf. 
eq. (Q). 



Ion 


z 


a [A] 


b [A] 


c [AVkJ/mol] 


Li+ 


0.87 


1.0155 


0.07321 


22.24 


Si 4 + 


2.40 


0.8688 


0.03285 


47.43 




-1.38 


2.0474 


0.17566 


143.98 


/o = 4.186 kjA^VoU 1 r c = 1.3 A 



ion sites and the diffusion path by comparing the results 
with those obtained in Sec. IIIII Finally, we give in Sec. I VII 
a summary of the results and an outlook to further re- 
search. 



II. MOLECULAR DYNAMICS SIMULATIONS 

We perform MD simulations of a Li2SiC>3 glass in the 
NVE ensemble at a temperature of 700 K with peri- 
odic boundary conditions. Newton's equations are solved 
using the Verlet algorithm with time step At — 1 fs. 
The computational domain is a cubic box of side length 
L = 16.68 A filled with 144 lithium, 72 silicon and 216 
oxygen ions. The box size was determined by performing 
simulations in the NPT ensemble at atmospheric pres- 
sure. It corresponds to material densities that match the 
experimental ones within 5% of error. 

Pair potentials of Gilbert-Ida typ o 25 i 26 describe the in- 
teractions between the species: 



where the parameters listed in table [I] have been 
optimized 23 and shown to give good agreement with ex- 
perimental datai 23 i 27 i 28 i 29 The first term in eq. Q is the 
Coulomb interaction with effective charge numbers Zi for 
Li, Si and O. The long-range Coulomb interaction with 
the image charges in the periodically continued copies 
of the simulation box is taken into account by stan- 
dard Ewald summation. A Born-Meyer type potential 
Ay exp(— r/Ay ) in the second term of eq. Q takes into 
account the repulsive short-range interactions, where the 
parameters et^ and bj appearing in JQl decompose Ay and 
Ay into values assigned to the interacting species. The 
last term in eq. Q is a dispersive interaction and present 
only for interactions between oxygen ions with distance 
larger than r c = 1.3 A^ 

We cooled down the system in the NVT ensemble us- 
ing a velocity scaling to reach the final temperature of 
T = 700 K and subsequently equilibrated the system in 
the NVE ensemble. For the analysis of the results we 
performed a simulation run over 40 ns. Figure^, shows 
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the pair correlation functions for Li-Li, Li-Si and Li-0 
that allow us to check if the minimal distances between 
pairs of ions used in the BV analysis (2.48 A for Li-Si 
and 1.7 A for Li-O, see below) fit to the potential model 
(|T|). The time-dependent mean square displacement of 
Li, O and Si ions is displayed in Figure The mean 
square displacement of the Li ions becomes normal dif- 
fusive at time scales larger than several hundred picosec- 
onds, while that of the Si and O ions practically stays 
constant over the whole time interval. On the time scale 
of 2 ns, for which most of the calculations for the de- 
termination of sites and diffusion paths are carried out, 
aging effects of the network structure can be neglected. 



III. IDENTIFICATION OF ION SITES AND 
THE DIFFUSION PATH 

During their motion the Li ions explore only parts of 
the host network. To identify the regions encountered 
by the Li ions and their favorable sites, the local num- 
ber density p(x) of Li ions is calculated. Dependent on a 
threshold value p t h, we define a path V(pth) — { x |p( x ) > 
Pth} as the region with p(x) > p t h, and perform a sub- 
sequent cluster analysis of the paths in dependence of 
Pth- To determine p(x) the simulation box is subdivided 
into a grid of cubic cells with spacing A and the time 
ti is registered, where it is occupied by a Li ion. The 
average density pi ~ p(x.i) in cell i is then calculated 
as pi — (ti/t s i m )A -3 , where t slm is the total simulation 
time. 

In order to identify regions of high probability of oc- 
cupation we determine clusters of connected (face shar- 
ing) cells i with pi > p t h by the Hoshen-Kopelman 
algorithm^ Figure [3 shows the number of clusters 
Aci(pth) in dependence of pth for three grid spacings A. 
All three curves have the same shape: Starting from large 
Pth, ^Vci(pth) increases with decreasing p t h, since an in- 
creasing number of local maxima in p(x) is identified. 
The strong dependence on the grid spacing at large p t h 
shows that the local maxima in p(x) are sharp. For larger 
A, less local maxima and consequently less clusters are 
found, since the average density pt in a cell containing a 
sharp local maximum of p(x) becomes smaller and can 
fall below p t h- As long as A is much smaller than the 
typical distance between the local maxima one expects 
that with decreasing p t h eventually all local maxima are 
resolved. This is indeed the case, as can be seen from the 
fact that the curves for different A pass the same maxi- 
mum at j\r™ ax ~ 165. By further decreasing p t h different 
clusters merge and A c i(p t h) decreases. 

The discussion makes clear that there exists a plateau 
N ci {pth) = N^ ax in the range p min < p t h < Pmax, which 
depends on the grid spacing (while the value N^ ax is not 
affected for sufficiently small A). Above p max some of 
the local maxima in p(x) are not resolved and below p m i n 
clusters coalesce (which does not exclude that some of the 
jV™ ax clusters contain more than one local maximum of 



p(x)). To identify clusters with possible Li sites one could 
choose any value of pth in the threshold range p mm < 
Pth < Pmax- This would not change the identity of the 
possible sites but their size. In order to cover as much 
volume as possible we use p* = p m i n ~ 0.56 A -3 for a 
unique definition of the clusters, which are candidates 
for sites (strictly speaking, one should use the value p m i n 
in the limit A — > 0, since p m in depends weakly on A). 
The subsequent analysis is carried out for the smallest 
spacing A = 0.139. 

Next we study properties of the clusters with respect 
to their assignment to sites. To this end we determine 
the volume V c \ = n ce n s A 3 of the clusters and the mean 
number density p c \ = X^ieci Pi/ n cc\\s of Li ions on them, 
and number the clusters in order of increasing p c \. n cc \\ s 
denotes the number of cells in the corresponding clus- 
ter. Figures |3^,b,c show p c \, V c \, and the occupation 
probability t cl /t sim = X«ed W*sim = p c i«ccii s A 3 for the 
numbered clusters. The nearly constant increase of the 
density is interrupted by a jump after cluster number 17. 
Also the other two quantities show a strong increase in 
this region. The first clusters are very small and have 
very low occupation times. After the jump, the volume 
and the occupation time lie on a main branch between 
0.6 and 0.35 A 3 for the volume and between 50 and 90% 
of the ratio of occupation time to simulation time. In 
between there are some outliers with less volume or less 
occupation time, but above the values of the first 17 clus- 
ters. With these results one faces two problems: (^Sum- 
ming up the occupation probability of all clusters, one 
finds that the total occupation probability i c iAsim of 
all clusters is 70% only, (ii) According to the behavior of 
all three shown quantities — very small occupation prob- 
ability, volume and occupation time — the first clusters 
before the jump do not fit to the physical picture which 
one has of a site. 

Problem (i) can be solved as follows. The criterion 
p(x) > p+ cuts off the outer fringe of a site, which is 
visited only a short time by a Li ion, but dynamically 
can be assigned to the site, since a Li ion in general re- 
turns to the cluster after entering its fringe. The clusters 
thus can be viewed to form the core of the sites. Prob- 
lem (ii) is more subtle. As suggested in ref. 0, one can 
require that a cluster should only be assigned to a site 
if the mean occupation of Li ions on it exceeds a cer- 
tain threshold. However, there can exist sites, which are 
visited only rarely but which still conform to the require- 
ment of a sufficiently large ratio between the residence 
time of a Li ion on the site and the time for its transition 
to a neighboring site. 

We therefore use a new criterion for the assignment of 
the clusters to sites by associating two times with each 
cluster, the total residence time i rcs and the total hopping 
time ihop to any of its neighbors. The hopping time is 
calculated as follows: When a Li ion leaves a cluster and 
enters a new one without returning, it performs a transi- 
tion between two clusters. The hopping time then is the 
time interval between the departure from the initial clus- 
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ter and the arrival at the target cluster. Summing over 
all transitions starting from a given cluster, one obtains 
the total hopping time from this cluster. The time inter- 
val between the entrance of a Li ion into a cluster and 
the onset of its transition to another cluster yields the 
residence time, and by taking the sum over all events we 
obtain t rcs £& Note that according to this definition the 
residence time includes events, where a Li ion leaves the 
cluster but returns to it before entering another cluster^. 

Figure 0Jl shows t res and ihop for the numbered clus- 
ters. While ihop fluctuates from cluster to cluster but 
never exceeds one tenth (0.2 ns) of the simulation time 
(2 ns) , i ros shows a more smooth variation with the clus- 
ter number as a consequence of a strong correlation to 
p c \. In particular, clusters with the smallest p c \ corre- 
spond to the one with lowest i res . The total residence 
probability t Ies /t S i m of the Li ions on the clusters is 
now 98%, which confirms that the problem with the low 
total occupation probability X^ci^ci/^sim is resolved by 
taking into account the fringes of the clusters. 

In order for a cluster to be assigned to a site, we now 
use the criterion that its residence time t s i m should be at 
least one order of magnitude larger than its total hopping 
time tho P r^ 



^re 
tho 



> 10. 



(2) 



Figure 0Jd shows £ re sAhop for the numbered clusters. 
Most clusters fulfill condition J5J, which demonstrates 
that the Li motion can be represented by a hopping dy- 
namics on a coarse grained time scale larger than ihop- 
Only 17 fail to satisfy condition (01, and almost all of 
those have small numbers corresponding to low p c \. This 
shows that there is hardly any difference here between a 
criterion based on a threshold number density and 
Accordingly, the sites are determined essentially by the 
equilibrium quantity p(x) (in the metastable basin after 
the cooling when disregarding slow non-equilibrium ag- 
ing processes of the network structure). In total, we find 
148 sites, which are only 2.8% more than the number of 
Li ions. This fraction of empty sites is slightly lower than 
that found in earlier studiesS*^ and supports the picture 
of an ion hopping with a small number of vacancies. Dif- 
ferent from the usual situation in crystalline systems in 
thermal equilibrium, the concentration of vacancies here 
should be considered as a result of the freezing process 
and will not change significantly with temperature below 
the glass transition. 

After having identified the sites, we determine the dif- 
fusion path. To this end we evaluate the percolation 
threshold^ p pe rc and the corresponding subset V pcrc = 
^(Ppcrc) of cells, on which the Li ions can diffuse across 
the system. We find p pcic = 0.030 A 3 , where P(p pC ic) 
covers 7.0% of the volume of the system (cf. table IIII|I . 
Confining the Li motion to exactly the critical perco- 
lation path at p = Ppcrc would yield anomalous sub- 
diffusion. However, in our case the percolation threshold 
is not sharp due to the finite system size and the number 



density of the Li ions takes into account only the mean 
occupation probabilities of the ions but not their thermal 
fluctuations. Even more important, these issues are not 
of particular relevance here, since we are interested in the 
spatial overlap with the BV paths (see Sec. II VII . With 
respect to this overlap the criticality plays no decisive 
role (since connectivity properties are not important for 
the overlap). 



IV. BOND VALENCE ANALYSIS 

The development of the bond valence approach has its 
origin in the search of correlations between bond lengths, 
chemical valence and binding energies for the chemical 
bondjSi and has been widely applied to crystals with co- 
valent and ionic bonds. The possibility to approximately 
describe both types of bonding using the same formalism 
makes it particularly useful for inorganic compounds with 
partially coval ent bo nds. A review of the BV method is 
given in refs. 20 21J. We describe the method here in 
connection with our application, which is the determi- 
nation of sites and diffusion paths for the Li ions based 
on the knowledge of the Si-0 network structure. For a 
Li ion to be accommodated in some region, its valence 
should be close to its "natural values", and additional 
constraints for its coordination number and minimal dis- 
tance to neighboring ions are to be satisfied (see below). 

According to Paulingip& each bond in a structure in- 
duces, due to its polarity, bond valences of opposite sign 
(partial charges in units of the elementary charge e) at 
the two atoms that it connects. In an equilibrated struc- 
ture, the bond valences of an ion induced by neighboring 
counterions should add to its ideal value (i.e. +1 for a Li 
ion) . Since the effective overlap of electronic orbitals typ- 
ically decreases exponentially with distance of the atomic 
nuclei, the partial valence Sj(x) of a Li ion at position x 
induced by a oxygen ion j at position x^ is determined 
by^ 



(x) = exp 



r - |x-Xj| 



(3) 



where ro is the ideal bond length and £ is the so-called 
softness parameter that determines how fast the bond 
valence varies with distance. The total valence V^(x) of a 
Li ion is computed as the sum of the bond valences Sj (x) : 



£ 



(4) 



Because of the exponential decay with distance the sum 
can be extended over all oxygen ions in the computa- 
tional domain (taking into account the minimum image 
convention for periodic boundary conditions). The de- 
composition J2J m t° bond valences allows one also to 
associate a coordination number C (x) of a Li number at 
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TABLE II: Parameters used in the BV analysis. 



£ [A] 


ro [A] 


Smin 


Cmin 


Cm ax 


RuO [A] 


i?LiSi [A] 


V 


0.516 


1.17096 


0.04 


4 


6 


1.70 


2.48 


0.814 



position x. This is denned by the number of oxygen ions 
contributing Sj(x) exceeding a threshold s min , 

C(x) = ^0( Sj -(x)- Smin ) , (5) 
j 

where #(.) is the Heaviside jump function = 1 for 

x > and zero else) 3 Here we choose s m i n as suggested 
by Brown^S The BV parameters ro, and £ are derived 
from a variety of crystalline phases^2i2S and are given in 
table |n| 

Figure shows the histogram of valences of the 144 
Li ions for an instantaneous equilibrated configuration. 
While in a crystalline structure the Li ions i at their 
equilibrium positions x^ have valences Vi = V(xj) close 
to the ideal value V;d = 1, we find a mean value VJnst. = 
0.821 significantly smaller than Vid in the (instantaneous) 
amorphous glass structure. An apparent under-bonding 
of the Li ions with average BV sums of ca. 0.9 has been 
found also in the BV analysis of RMC models obtained 
from scattering datai& and can be traced back to the 
non-crystalline structure of the glass. Since in the MD 
simulation the cooling rates are much larger than in ex- 
periments, the deviation from the ideal bonding situation 
in a crystal at thermal equilibrium is even larger. 

With respect to the application of the BV analysis 
to structural models obtained from RMC simulations, 
one should take into account that the information pro- 
vided by x-ray and neutron diffraction data corresponds 
to time-averaged positions of the network forming ions. 
This is due to the fact that the time for obtaining an 
evaluable signal is orders of magnitudes larger than a 
typical vibrational time. This means that a (non-unique) 
network configuration obtained in RMC modeling is a 
representative for a time-averaged density. We thus de- 
termined the mean positions of oxygen and silicon ions 
in a time interval of 2 ns and calculated the valences of 
Li ions with respect to these mean positions. For a Li 
ion placed at the centers Xj of the cells i described in 
the previous section ITTT1 fspacing A = 0.139 A), the va- 
lences Vi — V(xj) were calculated. Taking into account 
the probability of occupation p(xi)A 3 of the cells, we 
then determined the probability density p(V) of valences 
shown in Fig. |SJd. Its mean V = 0.814 is even slightly 
smaller than that for the instantaneous configuration, 
confirming that the ideal value Vld = 1 is not the pre- 
ferred one in the MD simulation of an amorphous glass 
structure. A small asymmetry is seen in p(V) with a 
steeper decrease from its maximum on the side of large 
valences V. This reflects the asymmetry of the two-body 
interaction potentials of Li with O ions, where small dis- 
tances corresponding to the repulsive part yield larger 



valences according to eq. Q, while large distances cor- 
responding to the attractive part yield smaller valences. 
Since the asymmetry is not pronounced, the mean value 
V is close to the value where p(V) attains its maximum. 

In the following we analyze the data with respect to the 
time-averaged network structure and consider V as the 
"optimal value" of the bond valence sum. The valence 
mismatch of a Li ion at position x is then given by 

d(x) = |V(x) - V\ . (6) 

The mean valence mismatch 12% calculated from p(V) 
in Fig. 03 is significantly larger than in crystalline struc- 
tures (where typical values are less than 5%). 

Following the BV method, a position x is accessible for 
a Li ion if the following conditions are fulfilled: 

(i) The distances of the Li cation to silicon cations 
must exceed a minimum distance 

min{|x-xf|}>iiLiSi, (7) 

j 

-^LiSi = 2.48 A is chosen to be equal to (or a lit- 
tle less than) the sum of the radii of the two ions 
in agreement with the correlation hole of the pair 
correlation shown in Fig. 

(ii) According to the "equal valence rule"— among the 
conceivable states with matching BV sum, the one 
with more symmetric bonds is energetically prefer- 
able. The simplest way of excluding that a match- 
ing BV sum is achieved by an unphysical strong 
asymmetric coordination shell is to define a mini- 
mum acceptable bond distance, 

min{|x-x°|} >i? Li0 . (8) 

3 

A value of -Ruo = 1.7 A is in agreement with the 
correlation hole in Fig. QJl and restricts a partial 
Li-0 bond valence to values smaller than 0.36. 

(Hi) The coordination number has to lie between a min- 
imal and maximal value, 

Cmin < C(x) < C max , (9) 

where C m i n = 4 and C max = 6. This condition has 
previously been assumed for sites only— 

(iv) The valence mismatch must be smaller than a 
threshold value dth, 

d(x) = \V(x) - V\ < dth ■ (10) 

As an alternative to the three conditions (ii-iv) one 
can use just one condition (v), which is a modification of 
(iv). It takes into account that a Li ion is better accom- 
modated to the local network environment at position x 
if the valences Sj (x) are more symmetrically distributed 
among the neighboring oxygen ions. The effect can be 
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TABLE III: Volume fractions (%) of the various subsets spec- 
ifying sites and diffusions paths. 



^sitcs 


^pcrc 


•psv 


-pBV 
' cn 


' cx+cn 


1.4 


7.0 


10.0 


83.0 


8.3 



described by defining an "ideal partial valence" Sid by 
CminSid = Vid i 18 ^ 4 corresponding to a Li ion that, when 
symmetrically connected to C m - ln oxygen ions, has the 
ideal total valence VJd- The deviation from the symmet- 
ric situation is quantified by the penalty function 




where the sum is taken over all oxygen ions and we choose 
/i = 3. A modified valence then is defined by d'(x) — 
d(x) + p(x) and we require: 

(v) The modified valence mismatch must be smaller 
than a threshold value dth, 

d , (x) = \V(x.)-V\+p(x)<d th . (12) 

When applying eq. I|12|) instead of eq. (|10fl conditions 
(H) and (Hi) are no longer needed. As a consequence, 
instead of the four parameters Ruo, Cmin, Cmax, and 
Smin the parameter \i enters the calculation. 

To perform the BV analysis, the computational do- 
main is divided into a grid of cells with spacing A anal- 
ogous to the procedure used in sec. Hill (A = 0.139 A). A 
cell is accessible if a Li ion placed at its center fulfills the 
requirements (i-iv), or (i,v) in the modified BV method^ 
The union of such cells is the BV path "Pevt^th) (for con- 
ditions (i-iv)) or the BV path V^ v (dth) 



V. COMPARISON OF ION SITES AND 
DIFFUSION PATHS 

We first clarify the relevance of the purely geometric 
constraints (i,ii) in combination with the coordination 
number condition (Hi), without taking into account the 
bond valence condition (iv). Accordingly we distinguish 
between the subsets defined by applying the constraints 
separately: is the subset given by the geometric 

exclusions (i,H), while refers to the condition im- 

posed solely on the coordination number (Hi). The sub- 
set defined by the combined conditions (i-iii) is denoted 

The geometric conditions (i,ii) already limit the frac- 
tion of available space for the mobile ions to 10% for 
the time-averaged network structure, see table IIIII By 
contrast, constraint (Hi) alone for s m j n = 0.04 is a weak 
condition, which would allow a Li ion to occupy 83% of 



TABLE IV: Comparison of the subsets determined by the 
B V analysis according to criteria ( i-iii ) with the sites and dif- 
fusion paths identified from the MD simulations. Upper part: 
Sensitivities 7/)('Pj 3V |'P t ); lower part: specificities ^>('P*|'P;P V ). 
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0.12 


0.02 


0.13 


T^pcrc 


0.43 


0.08 


0.48 



the space. The combined conditions (i-iii) reduce the ac- 
cessible space to 8.3%, corresponding to an uncorrelated 
behavior. If conditions (ii-iv) are replaced by condition 
(v) in the limit d t h — > oo the available space is restricted 
to 21.3% of the volume of the system (see also Fig. 
One should expect that this accessible space entails the 
ionic sites and the diffusion path identified in sec. IIIII 

To quantify the quality of agreement between the sub- 
sets P BV (P BV = 7> BV , P BV , or P^ +cn ) with the sub- 
sets V* obtained from the MD simulations (P+ = V s ites 
or "Ppcrc), we define two quantities, the sensitivity and 
specificity. The sensitivity is the conditional probability 
- 0('P BV |'P+) that a cell belonging to one of the subsets 
identified in the MD simulations belongs to one of the 
subsets of the BV analysis. Conversely, the specificity of 
the BV analysis is quantified by calculating the condi- 
tional probabilities VCP*|P BV ) — 

Table IIVI summarizes the results. Let us in particu- 
lar consider the case where conditions (i-iii) are applied 
(corresponding to P^+ cn )- A high sensitivity of 79% is 
reached for the sites. However, from table lllll we see that 
the volume fraction of V^+cn (8,3%) is by a factor of 
about 6 larger than that of Twites (1,4%). Accordingly, 
the specificity ^i(P S \tes\P^ +cn ) = 13% is rather low. The 
sensitivity for the diffusion path is 57% (Ppcrc)- It is 
significantly lower than that for the sites, since the path 
contains regions with very low occupation probability (cf. 
Sec. IIII|I , which often violate the geometric constraints 
(i,ii). 

The question is, whether the specificity can be im- 
proved without significant reduction in the sensitivity by 
applying condition (iv) in addition to (i-iii). The volume 
fraction v^v of the BV path as a function of the thresh- 
old mismatch dth is shown in Fig. (solid line). For small 
dth, v bv increases linearly and for dth ^ 0.3 approaches 
the value 8.3% imposed by the constraints (i-iii). Fig- 
ures [7^ and b show how the sensitivity and specificity 
vary with dth with respect to the sites and the diffusion 
path, respectively. The sensitivities in Figs. [7^,b start 
to saturate for dth ^ 0.2 — 0.3, a value in fair agree- 
ment with typical choices used in BV analyses of RMC 
models. 42 Close to dth — 0.1, vbv — 5% is by about 40% 
smaller than the saturation value 8.3% for e? t h —> oo (see 
Fig. |SJ. However, this reduction is not strong enough to 
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yield a significant improvement of the specificities. Even 
if one would take a large loss in the sensitivity by choos- 
ing dth very small, the gain in the specificity remains low. 
In summary we find that the inclusion of criterion (iv) 
does not yield a substantial improvement and that the 
essential part of the obtained agreement is due to the 
geometric constraints 

To find an optimal value for dth with respect to both 
sensitivity and specificity, we use the maximum of Co- 
hen's kappa valued. The kappa value for the sets A = 
P BV (d) and B = is defined by 

= ( P (A nB)+p(AnB))- ( P (A)p(B) + p(A)p(B)) 
l-(p(A) P (B)+p(A)p(B)) 

(13) 

where p(.) denote the probabilities of the corresponding 
sets (for example, p(A (~l B) is the probability that a cell 
belongs to both the sets A and B). A value of k = 1 
denotes complete agreement between the two sets, while 
k = corresponds to a random overlap of the two sets 
(negative values indicate an anti-correlation). For the 
sites the maximum of k occurs at dth = 0.09 and is only 
25%. This is caused by the low specificity, which varies 
only slowly with dth. At dth = 0.09 we find a specificity 
of 17%, while the sensitivity is 59%. For the path the 
maximum of k occurs at dth = 0.23 and is 48%. The 
sensitivity and specificity attain values of 54% and 50% 
at this maximum. In summary, the quality of agreement 
for the path is promising, while that for the sites is not 
yet satisfactory. 

The sensitivity for the diffusion path can however be 
improved and the influence of the BV sum mismatch in- 
creased by using condition (v) instead of conditions (ii- 
iv), i.e. by considering the path T-Bv(^th)- For this path 
in comparison with Vperc we plot in Fig.[S]the sensitivity 
and the specificity in dependence of <ith (see eq. I|12[l ■ For 
large dth, the sensitivity reaches values up to 90%, which 
are significantly higher than the 57% in Fig. 0d. The 
specificity for small dth has values comparable to that in 
Fig. 03, and then decreases to smaller values for larger 
dth- As a function of d t h, k reaches a maximum of 48% 
at dth = 0.22. Further enhancements may be achieved 
by slight adjustments of the minimum distance i?LiSi or 
by replacing criterion (i) by a penalty function as well, 
which would mean that the oversimplifying hard sphere 
exclusion radius criteria may adversely affect the 

achievable level of agreement. 

Finally, we test the potential of the BV method to 
identify sites (despite the low specificity). To this end we 
performed a cluster analysis with respect to d t h, analo- 
gous to the one carried out in sec. IIIII with respect to p c \. 
Different from the behavior found in Fig. [21 the number 
N® v of BV clusters shown in Figgis very large already 
for small dth and then decreases rapidly due to coales- 
cence of clusters. Accessible cells start to percolate at 
d c = 0.09. The reason for this behavior is that the va- 
lence mismatch d[x) is a rather rapidly varying function. 
Already for small d t h there exists a large number of small 



disconnected clusters that soon merge together to form 
a percolating cluster. While the majority of vbv belongs 
to this percolation cluster the fluctuations in d(x) pre- 
vent N BV to assume comparatively small values for large 
dth- Several hundred clusters are found for dth — 0.1. 
Almost all of these are very small: 80% consist of only a 
single cell and only 8 clusters contain more than 100 cells 
(i.e. are comparable to the size of sites). The percolation 
cluster includes 95% of «Bv(^th = 0.1). As a consequence 
the cluster analysis is not successful for identifying ionic 
sites with respect to d t h~ 

It is interesting, however, that the BV analysis can 
provide a good estimate of the number of sites (not their 
location in space). It turns out to be useful to include 
the coordination number constraint (Hi) in addition to 
the conditions (i,v) in this case. For varying mismatch 
threshold dth the cells belonging to the corresponding 
path "P B v(dth) are filled row by row under the condition 
that two centers have a distance larger than the minimal 
Li-Li distance i?LiLi = 2.62. The number of sites found in 
this way is shown in Fig. Ellas a function of the mismatch 
threshold dth (lower panel). For large dth, it approaches 
the number of clusters found by the Hoshen/Kopelman 
analysis in Sect. lIIll On the other hand, we can take the 
estimated number of sites at the optimal bond valence 
mismatch dth = 0.22 given by the maximum of Cohen's 
k value with respect to the sets "P s ites arid V^y(dth)- The 
result displayed in Fig. |S| yields 151 sites. This value is 
surprisingly close to the 148 sites found in Sec. IIIII but 
studies of further systems are required to test whether 
such a prediction is generally possible. 

In view of the success of the BV method to estimate 
the number of sites, we undertook further attempts to 
improve the level of agreement for the localization of the 
Li sites. If one distributes Li ions more randomly on 
the path PevO^th) with the same minimal distance con- 
straint i?LiLi = 2.62 A^ instead of doing it row by row as 
described above, the positioning of centers of sites is no 
longer biased to the rim of the BV path. In this way, we 
obtain 165 BV clusters for dth = 0.14. Considering now 
as BV sites the spheres with ion radius Ruhi/2 — 1.31 A, 
we found that 75% of these sites have at least one cell in 
common with one of the MD sites found in Sec. EH Still 
this is not enough for a reliable localization of Li sites. 



VI. SUMMARY AND PERSPECTIVES 

We have identified sites and diffusion paths for the 
mobile Li ions in molecular dynamics simulations of a 
Li2Si03 glass. The identification was based on a cluster 
analysis of regions with high Li number density p(x). For 
the clusters to be assigned to sites, we chose the condition 
that the total residence time of a Li ion on it (including 
intermediate escape to fringe regions) is ten times larger 
than the total hopping time to any of its neighboring 
clusters. Using this criterion, a very low concentration 
2,8% of vacant sites was found. An attempt to identify 
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the diffusion path and sites by bond valence sums calcu- 
lations was carried out. The comparison showed that the 
core of sites and parts of the diffusion path (up to about 
50%) are captured. Conversely, the BV method was not 
suitable to distinguish between regions of high and low 
p(x), and accordingly not suitable to determine the sites. 

When dealing with the MD reference data of this work 
it proves to be necessary to adapt the BV parameters 
derived from experimental diffraction data to the MD 
force model. We showed that it is advantageous to opti- 
mize or replace the oversimplifying hard sphere exclusion 
radii for application of the BV analysis to the MD model. 
Moreover, an improvement of the BV analysis seems to 
be possible by adjusting the BV parameters ro and £ to 
the MD force model instead of only adjusting Vm to V. 

In searching for a powerful method to relate struc- 
tural properties of the disordered network structure to 
transport properties of the mobile ions a further ap- 
proach could be based on an effective potential U c g (x) oc 
— fcBTlnp(x). Taking this effective potential one could 
perform a critical path analysis to determine the acti- 
vation energy for the long-range ion mobility (see e.g. 
ref . l46f) . This approach would incorporate Coulomb in- 
teraction effects between the mobile ions in a mean field 
type approximation. Preliminary studies by us show that 
such Coulomb effects are important for the formation of 
the sites. With respect to the BV method we have also 
performed preliminary studies to evaluate the degree of 
correlation between C/ e ff(x) and the BV mismatch <i(x) 
(or d'(x)). It turned out that BV method cannot clearly 
distinguish between regions of high and low p(x) and ac- 
cordingly it is not suitable to locate the sites. Besides 



the described shortcomings of a particular choice of pa- 
rameter values this is mainly caused by neglecting the 
long-range Li-Li Coulomb repulsions. 

The precise decomposition of the ion trajectories into 
residence and transition parts described in Sec. lIIII morc- 
over may allow one to bridge the time scale gap be- 
tween molecular dynamics and coarse-grained Monte- 
Carlo simulations. Given two neighboring sites, many 
transitions of mobile ions can be followed and the av- 
erage transition rate calculated. The elementary rates 
for the possible transitions can then be used in a subse- 
quent Monte-Carlo simulation. However, such procedure 
is not particularly useful if one is not able to take into 
account the temperature dependence of the elementary 
rates. One possibility is to use MD data at various high 
temperatures, and to try to calculate from them activa- 
tion energies for the elementary rates (if these follow an 
Arrhenius type behavior). With less effort, one could use 
again the density p(x) and the effective potential derived 
from it to determine the energetics of a coarse-grained 
hopping model. Further investigations in this direction 
have to be undertaken in order to establish a working 
multi-scale modeling of these complex amorphous sys- 
tems. 
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FIG. 1: (a) Pair correlation function of Li-Li, Li-Si and Li-O 
and (b) time-dependent mean square displacements of Li, O 
and Si ions at T = 700 K. 
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FIG. 2: Number of clusters N c \ in dependence of p t h for 
three different grid spacings A. Note that the mean number 
density of Li ions is 144/(16.68A) 3 ~ 0.03lA" 3 only, which 
in comparison with the relevant scale of p t h implies that the 
Li ions are strongly localized on the clusters. 
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FIG. 3: (a) Mean number density p c \ of Li ions on the 165 
clusters identified by the Hoshen-Kopelman algorithm; (b) 
Volume Vd and (c) Li occupation probability t c i/t sim of the 
clusters. The clusters are sorted according to increasing p c \. 
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FIG. 4: (a) Total residence and hopping time normalized 
with respect to the simulation time t s im, and (b) their ratio 
for the 165 clusters determined by the Hoshen-Kopelman al- 
gorithm. The dashed line marks the threshold above which 
clusters are identified as sites (cf. eq.|2Jl. 
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FIG. 5: (a) Histogram of of bond valence sums for an instan- 
taneous configuration of the equilibrated Li2SiOa glass; (b) 
Probability density of bond valence sums with respect to the 
time-averaged structure of network forming O and Si ions. 




FIG. 6: Volume fraction vbv of the BV path as a functions 
of the mismatch threshold d th . 
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FIG. 7: Quality of agreement between (a) the subsets "Psites 
and 7>Bv(cJth), and (b) the subsets "P pa th and 7 3 Bv(dth)- 




FIG. 8: Quality of agreement between the percolation path 
Pporc and Psvidth) as a function of the BV mismatch thresh- 
old d t h- 
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FIG. 9: Number of BV clusters as a function of the threshold 
mismatch dth- 




FIG. 10: Lower panel: Number of possible site centers on 
the BV path T-Bv(^th) as determined by the "row-by-row al- 
gorithm" described in the text; Upper panel: Cohen's kappa 
value for the path 7- > Bv(^th) in comparison with Ppcrc (re- 
drawn from Fig. [SJ- The number of estimated sites (151) at 
the optimal « is indicated by the solid lines. The dashed line 
marks the number of Li ions (144) and the dotted line the 
number of sites (148) found in Sec. IIIII 



